library(dplyr); library(ggplot2); library(ggpattern)

for(plot in "Figure1-a"){
  ggdata <- read.csv("Data-Fig1(a).csv")
  ggdata$River <- rep(c("c","b","a"),27) #c=East , b=South, a=West
  ggdata$Year <- as.character(ggdata$Year) %>% as.integer()
  
  ggplot(data=ggdata, aes(x=Year, y=number_of_station, 
                          pattern=River, pattern_angle=River, pattern_spacing=River))+
    geom_col_pattern(fill="white", color="black", width=0.6, linewidth=0.5, position="stack",
                     pattern_density=0.3, pattern_fill="black", pattern_color="black")+
    geom_vline(xintercept=2011.5, linetype="dashed", color="red2", linewidth=1)+
    
    scale_pattern_discrete(choices = c("stripe", "circle", "none"), labels=c("West", "South", "East"))+
    scale_pattern_angle_discrete(range = c(45,0,0), labels=c("West", "South", "East"))+
    scale_pattern_spacing_discrete(range = c(0.03, 0.01, 0), labels=c("West", "South", "East"))+
    
    scale_x_continuous(breaks = seq(1990, 2021, 5), expand=c(0.01,0.01))+
    scale_y_continuous(breaks = seq(0, 80, 5), expand=c(0.01,0.01))+
    coord_cartesian(xlim = c(1988.5, 2021.5), ylim = c(0, 22.5)) + 
    theme_bw()+ xlab(NULL)+ ylab("")+ 
    ggtitle("(a) River CQ Paired Stations") +
    theme(panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(),
          panel.grid.major.y = element_blank(), panel.grid.minor.y = element_blank(),
          axis.text.x = element_text(margin = margin(t=5), face = 1, size=18, color="Black"),
          axis.text.y.left = element_text(margin=margin(l=5,r=3), size=20, color="black"),
          axis.title.y.left = element_text(margin = margin(r=3), size=20, color="black"), 
          axis.text.y.right = element_text(margin=margin(r=11), size=20, color="white"),
          axis.title.y.right  = element_text(margin = margin(l=0), size=20, color="white"), 
          plot.title = element_text(margin=margin(b=10), hjust = 0, size = 22),
          axis.ticks.length.y.left = unit(1.5,"mm"), axis.ticks.length.y.right = unit(0,"mm"), 
          axis.ticks.length.x = unit(1.5,"mm"),
          axis.ticks.y.left = element_line(linewidth = 1.2, color="black"), 
          axis.ticks.y.right = element_line(linewidth = 0, color="white"), 
          axis.ticks.x = element_line(linewidth = 1.2), 
          legend.position = "none",#legend.position = c(0.1,0.7),
          panel.border = element_rect(color="black", fill=NA, linewidth=1))
  
  #ggsave("Figure 1 - (a).png", units="in", width=12, height=4, dpi=1200)
}
for(plot in "Figure1-b"){
  ggdata <- read.csv("Data-Fig1(b)1.csv")
  ggdata1 <- read.csv("Data-Fig1(b)2.csv")

  ggdata1$Inshore <- c(rep("c",25),rep("b",25),rep("a",25)) #c=East , b=South, a=West
  ggdata1$Year <- as.character(ggdata1$Year) %>% as.integer()
  dual_axis_ratio = 3.33
  ggdata1$number_of_stations <- ggdata1$number_of_stations / dual_axis_ratio
  
  ggplot()+
    geom_point(data=ggdata, aes(x=Year, y=number_of_stations), fill="white", color="blue3", size=3, shape=21, stroke=2)+
    geom_col_pattern(data=ggdata1, aes(x=Year, y=number_of_stations, 
                                       pattern=Inshore, pattern_angle=Inshore, pattern_spacing=Inshore),
                     fill="white", color="black" , width=0.7, linewidth=0.5, position="stack",
                     pattern_density=0.3, pattern_fill="black", pattern_color="black")+
    
    scale_pattern_discrete("Mariginal Sea",
                           choices = c("stripe", "circle", "none"), labels=c("West", "South", "East"))+
    scale_pattern_angle_discrete("Mariginal Sea",
                                 range = c(45,0,0), labels=c("West", "South", "East"))+
    scale_pattern_spacing_discrete("Mariginal Sea",
                                   range = c(0.03, 0.01, 0), labels=c("West", "South", "East"))+
    
    scale_x_continuous(breaks = seq(1990, 2021, 5), expand=c(0.01,0.01))+
    scale_y_continuous(breaks = seq(0, 80, 20), expand=c(0.01,0.01),
                       sec.axis = sec_axis(~.*dual_axis_ratio, name="Marginal Sea", 
                                           breaks=seq(0, 250, 50)))+
    coord_cartesian(xlim = c(1988.5, 2021.5), ylim = c(0, 82))+
    theme_bw()+ xlab(NULL)+  ylab("River")+ 
    ggtitle("(b) Water Quality Monitoring Sites") +
    theme(panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(),
          panel.grid.major.y = element_blank(), panel.grid.minor.y = element_blank(),
          axis.text.x = element_text(margin = margin(t=5), face = 1, size=18, color="Black"),
          axis.text.y.left = element_text(margin=margin(l=5,r=3), size=20, color="blue3"),
          axis.title.y.left = element_text(margin = margin(r=3), size=20, color="blue3"), 
          axis.text.y.right = element_text(margin=margin(r=11), size=20, color="black"),
          axis.title.y.right  = element_text(margin = margin(l=0), size=20, color="black"), 
          plot.title = element_text(margin=margin(b=10), hjust = 0, size = 22),
          axis.ticks.length.y.left = unit(1.5,"mm"), axis.ticks.length.y.right = unit(1.5,"mm"), 
          axis.ticks.length.x = unit(1.5,"mm"),
          axis.ticks.y.left = element_line(linewidth = 1.2, color="blue3"), 
          axis.ticks.y.right = element_line(linewidth = 1.2, color="black"), 
          axis.ticks.x = element_line(linewidth = 1.2), 
          legend.position = "none",
          panel.border = element_rect(color="black", fill=NA, linewidth=1))
  
  #ggsave("Figure 1 - (b).png", units="in", width=12.36, height=4, dpi=1200)
}

for(plot in "Figure3-a"){
 #[Precipitation]
  ggdata <- read.csv("Data-Fig3.csv")
  ggplot()+
    geom_line(data=ggdata, aes(x=Year, y=Precipitation, linetype=Region), color= "black", linewidth=1)+
    geom_point(data=ggdata, aes(x=Year, y=Precipitation, shape=Region), 
               color= "black", fill="white", size=5, stroke=2)+
    
    scale_linetype_manual(values=c("solid", "solid", "solid", "solid"), 
                          breaks=c("Korea","West","South","East"))+
    scale_shape_manual(values = c(19,21,24,22), 
                       breaks=c("Korea","West","South","East"))+
    scale_x_continuous(breaks = seq(2012, 2021, 2))+
    scale_y_continuous(breaks = seq(0, 20, 5))+
    coord_cartesian(xlim = c(2012, 2021), ylim = c(0, 11.5)) + 
    theme_bw()+ xlab(NULL)+ ylab("[m/yr]")+
    ggtitle("(a) Precipitation (P)")+
    theme(panel.grid.major.x = element_line(colour="black", linewidth=0.1),
          panel.grid.minor.x = element_line(colour="black", linewidth=0.1),
          panel.grid.major.y = element_line(colour="black", linewidth=0.1),
          panel.grid.minor.y = element_line(colour="black", linewidth=0.1),
          axis.text.x = element_text(margin = margin(t=5), face = 1, size=25, color="black", angle=0),
          axis.text.y.left = element_text(margin=margin(l=4,r=5), size=25, color="black"),
          axis.title.y.left = element_text(margin = margin(r=5), size=25, color="black"), 
          plot.title = element_text(margin=margin(b=10), hjust = 0, size = 25),
          axis.ticks.length.y = unit(2,"mm"), axis.ticks.length.x = unit(2,"mm"),
          axis.ticks.y = element_line(linewidth = 1, color="black"), 
          axis.ticks.x = element_line(linewidth = 1, color="black"), 
          legend.position = c(0.275,0.9), legend.background = element_blank(),
          legend.title = element_blank(), legend.text = element_text(size = 25),
          legend.key = element_rect(color=NA, fill=NA, linewidth=1),
          legend.key.size=unit(10,"mm"),legend.direction = "horizontal",
          panel.border = element_rect(color="black", fill=NA, linewidth=1))
  
  #ggsave("Figure 3 - (a).png", units="in", width=11.59, height=3.6, dpi=1200)
}
for(plot in "Figure3-b"){
  #[Streamflow]
  ggdata <- read.csv("Data-Fig3.csv")
  ggplot()+
    geom_line(data=ggdata, aes(x=Year, y=Streamflow, linetype=Region), color= "black", linewidth=1)+
    geom_point(data=ggdata, aes(x=Year, y=Streamflow, shape=Region), color= "black", fill="white", size=5, stroke=2)+
    
    scale_linetype_manual(values=c("solid", "solid", "solid", "solid"), 
                          breaks=c("Korea","West","South","East"))+
    scale_shape_manual(values = c(19,21,24,22),  
                       breaks=c("Korea","West","South","East"))+
    scale_x_continuous(breaks = seq(2012, 2021, 2))+
    scale_y_continuous(breaks = seq(0, 20, 5))+
    coord_cartesian(xlim = c(2012, 2021), ylim = c(0, 13)) + 
    theme_bw()+ xlab(NULL)+ ylab("[Gm"^3~"/yr]")+
    ggtitle("(b) Streamflow (Q)")+
    theme(panel.grid.major.x = element_line(colour="black", linewidth=0.1),
          panel.grid.minor.x = element_line(colour="black", linewidth=0.1),
          panel.grid.major.y = element_line(colour="black", linewidth=0.1),
          panel.grid.minor.y = element_line(colour="black", linewidth=0.1),
          axis.text.x = element_text(margin = margin(t=5), face = 1, size=25, color="black", angle=0),
          axis.text.y.left = element_text(margin=margin(l=0,r=5), size=25, color="black"),
          axis.title.y.left = element_text(margin = margin(r=8), size=25, color="black"), 
          plot.title = element_text(margin=margin(b=10), hjust = 0, size = 25),
          axis.ticks.length.y = unit(2,"mm"), axis.ticks.length.x = unit(2,"mm"),
          axis.ticks.y = element_line(linewidth = 1, color="black"), 
          axis.ticks.x = element_line(linewidth = 1, color="black"), 
          legend.position = "none", legend.background = element_blank(),
          legend.title = element_blank(), legend.text = element_text(size = 25),
          legend.key = element_rect(color=NA, fill=NA, linewidth=1),
          legend.key.size=unit(10,"mm"),legend.direction = "horizontal",
          panel.border = element_rect(color="black", fill=NA, linewidth=1))
  
  #ggsave("Figure 3 - (b).png", units="in", width=12.1, height=3.6, dpi=1200)
}
for(plot in "Figure3-c"){ 
  #[TN Flux]  
  ggdata <- read.csv("Data-Fig3.csv")
  ggplot()+
    geom_line(data=ggdata, aes(x=Year, y=TN, linetype=Region), color= "black", linewidth=1)+
    geom_point(data=ggdata, aes(x=Year, y=TN, shape=Region), color= "black", fill="white", size=5, stroke=2)+
    
    scale_linetype_manual(values=c("solid", "solid", "solid", "solid"), 
                          breaks=c("Korea","West","South","East"))+
    scale_shape_manual(values = c(19,21,24,22), 
                       breaks=c("Korea","West","South","East"))+
    scale_x_continuous(breaks = seq(2012, 2021, 2))+
    scale_y_continuous(breaks = seq(0, 100, 20))+
    coord_cartesian(xlim = c(2012, 2021), ylim = c(0, 60)) + 
    theme_bw()+ xlab(NULL)+ ylab("[Gg/yr]")+
    ggtitle("(c) Total Nitrogen (TN) Flux")+
    theme(panel.grid.major.x = element_line(colour="black", linewidth=0.1),
          panel.grid.minor.x = element_line(colour="black", linewidth=0.1),
          panel.grid.major.y = element_line(colour="black", linewidth=0.1),
          panel.grid.minor.y = element_line(colour="black", linewidth=0.1),
          axis.text.x = element_text(margin = margin(t=5), face = 1, size=25, color="black", angle=0),
          axis.text.y.left = element_text(margin=margin(l=5,r=5), size=25, color="black"),
          axis.title.y.left = element_text(margin = margin(r=8), size=25, color="black"), 
          plot.title = element_text(margin=margin(b=10), hjust = 0, size = 25),
          axis.ticks.length.y = unit(2,"mm"), axis.ticks.length.x = unit(2,"mm"),
          axis.ticks.y = element_line(linewidth = 1, color="black"), 
          axis.ticks.x = element_line(linewidth = 1, color="black"), 
          legend.position = "none", legend.background = element_blank(),
          legend.title = element_blank(), legend.text = element_text(size = 25),
          legend.key = element_rect(color=NA, fill=NA, linewidth=1),
          legend.key.size=unit(10,"mm"),legend.direction = "horizontal",
          panel.border = element_rect(color="black", fill=NA, linewidth=1))
  
  #ggsave("Figure 3 - (c).png", units="in", width=12, height=3.6, dpi=900)
}
for(plot in "Figure3-d"){ 
  #[TP Flux] 
  ggdata <- read.csv("Data-Fig3.csv")
  ggplot()+
    geom_line(data=ggdata, aes(x=Year, y=TP, linetype=Region), color= "black", linewidth=1)+
    geom_point(data=ggdata, aes(x=Year, y=TP, shape=Region), color= "black", fill="white", size=5, stroke=2)+
    
    scale_linetype_manual(values=c("solid", "solid", "solid", "solid"), 
                          breaks=c("Korea","West","South","East"))+
    scale_shape_manual(values = c(19,21,24,22), 
                       breaks=c("Korea","West","South","East"))+
    scale_x_continuous(breaks = seq(2012, 2021, 2))+
    scale_y_continuous(breaks = seq(0, 10, 0.5))+
    coord_cartesian(xlim = c(2012, 2021), ylim = c(0, 1.8)) + 
    theme_bw()+ xlab(NULL)+ ylab("[Gg/yr]")+
    ggtitle("(d) Total Phosphorus (TP) Flux")+
    theme(panel.grid.major.x = element_line(colour="black", linewidth=0.1),
          panel.grid.minor.x = element_line(colour="black", linewidth=0.1),
          panel.grid.major.y = element_line(colour="black", linewidth=0.1),
          panel.grid.minor.y = element_line(colour="black", linewidth=0.1),
          axis.text.x = element_text(margin = margin(t=5), face = 1, size=25, color="black", angle=0),
          axis.text.y.left = element_text(margin=margin(l=0,r=4), size=25, color="black"),
          axis.title.y.left = element_text(margin = margin(r=8), size=25, color="black"), 
          plot.title = element_text(margin=margin(b=10), hjust = 0, size = 25),
          axis.ticks.length.y = unit(2,"mm"), axis.ticks.length.x = unit(2,"mm"),
          axis.ticks.y = element_line(linewidth = 1, color="black"), 
          axis.ticks.x = element_line(linewidth = 1, color="black"), 
          legend.position = "none", legend.background = element_blank(),
          legend.title = element_blank(), legend.text = element_text(size = 25),
          legend.key = element_rect(color=NA, fill=NA, linewidth=1),
          legend.key.size=unit(10,"mm"),legend.direction = "horizontal",
          panel.border = element_rect(color="black", fill=NA, linewidth=1))
  
  #ggsave("Figure 3 - (d).png", units="in", width=12, height=3.6, dpi=900)
}
for(plot in "Figure3-e"){ 
  #[Chl_a]
  ggdata <- read.csv("Data-Fig3.csv")
  ggplot()+
    geom_line(data=ggdata, aes(x=Year, y=Chl_a, linetype=Region), color= "black", linewidth=1)+
    geom_point(data=ggdata, aes(x=Year, y=Chl_a, shape=Region), color= "black", fill="white", size=5, stroke=2)+
    
    scale_linetype_manual(values=c("solid", "solid", "solid", "solid"), 
                          breaks=c("Korea","West","South","East"))+
    scale_shape_manual(values = c(19,21,24,22), 
                       breaks=c("Korea","West","South","East"))+
    scale_x_continuous(breaks = seq(2012, 2021, 2))+
    scale_y_continuous(breaks = seq(0, 10, 0.1))+
    coord_cartesian(xlim = c(2012, 2021), ylim = c(-0.01, 0.41)) + 
    theme_bw()+ xlab(NULL)+ ylab("[Gg/yr]")+
    ggtitle("(e) Chlorophyll-a (Chl-a) Flux")+
    theme(panel.grid.major.x = element_line(colour="black", linewidth=0.1),
          panel.grid.minor.x = element_line(colour="black", linewidth=0.1),
          panel.grid.major.y = element_line(colour="black", linewidth=0.1),
          panel.grid.minor.y = element_line(colour="black", linewidth=0.1),
          axis.text.x = element_text(margin = margin(t=5), face = 1, size=25, color="black", angle=0),
          axis.text.y.left = element_text(margin=margin(l=0,r=4), size=25, color="black"),
          axis.title.y.left = element_text(margin = margin(r=8), size=25, color="black"), 
          plot.title = element_text(margin=margin(b=10), hjust = 0, size = 25),
          axis.ticks.length.y = unit(2,"mm"), axis.ticks.length.x = unit(2,"mm"),
          axis.ticks.y = element_line(linewidth = 1, color="black"), 
          axis.ticks.x = element_line(linewidth = 1, color="black"), 
          legend.position = "none", legend.background = element_blank(),
          legend.title = element_blank(), legend.text = element_text(size = 25),
          legend.key = element_rect(color=NA, fill=NA, linewidth=1),
          legend.key.size=unit(10,"mm"),legend.direction = "horizontal",
          panel.border = element_rect(color="black", fill=NA, linewidth=1))
  
  #ggsave("Figure 3 - (e).png", units="in", width=12, height=3.6, dpi=900)
}

for(plot in "Figure4-a"){ 
  # [TN]
  # Korea
  bardata <- read.csv("Data-Fig4.csv") %>% filter(Species=="TN") %>% filter(Region=="Korean")
  ggplot()+
    geom_col_pattern(data=bardata, aes(x=Flux_Sum_percent, y=Type, pattern=Type, pattern_angle=Type),
                     fill="white",color="black", width=0.7, linewidth=2,
                     pattern_density=0.05, pattern_fill="Grey50", pattern_color="Grey50",
                     pattern_spacing=c(0.04,0.02,0.05, 1))+
    scale_y_discrete(limits=c("error", "dQdt_C", "dCdt_Q", "dCQdt"))+
    coord_cartesian(xlim = c(-100,50))+ 
    theme_void()+
    theme(legend.position = "none", axis.ticks.length = unit(0,"mm"),
          axis.text = element_blank(), axis.ticks = element_blank())+
    annotate("text", x = 34, y = bardata$Type, angle = -90, color= "black", fontface=2, size=10,
             label = c("F*", "C*", "Q*", "E*"))+
    annotate("text", x = 22, y = bardata$Type, angle = -90, color= "black", fontface=2, size=10,
             label = paste(round(bardata$Flux_Sum_percent,1),"%",sep=""))+
    annotate("text", x = 10, y = bardata$Type, angle = -90, color= "black", fontface=2, size=8,
             label = paste("(",round(bardata$Flux_Sum,0),")",sep=""), )+
    geom_vline(xintercept=0, color="Gray20", linewidth = 2)
  #ggsave("Figure4 - (a) Korea.png", units="in", width=6, height=6, dpi=1500)
}
for(plot in "Figure4-b"){
  #[TP]
  # Korea
  bardata=read.csv("Data-Fig4.csv") %>% filter(Species=="TP") %>% filter(Region=="Korean")
  ggplot()+
    geom_col_pattern(data=bardata, aes(x=Flux_Sum_percent, y=Type, pattern=Type, pattern_angle=Type),
                     fill="white",color="black", width=0.7, linewidth=2,
                     pattern_density=0.05, pattern_fill="Grey50", pattern_color="Grey50",
                     pattern_spacing=c(0.04,0.02,0.05, 1))+
    scale_y_discrete(limits=c("error", "dQdt_C", "dCdt_Q", "dCQdt"))+
    coord_cartesian(xlim = c(-100,50))+ 
    theme_void()+
    theme(legend.position = "none", axis.ticks.length = unit(0,"mm"),
          axis.text = element_blank(), axis.ticks = element_blank())+
    annotate("text", x = 34, y = bardata$Type, angle = -90, color= "black", fontface=2, size=10,
             label = c("F*", "C*", "Q*", "E*"))+
    annotate("text", x = 22, y = bardata$Type, angle = -90, color= "black", fontface=2, size=10,
             label = paste(round(bardata$Flux_Sum_percent,1),"%",sep=""))+
    annotate("text", x = 10, y = bardata$Type, angle = -90, color= "black", fontface=2, size=8,
             label = paste("(",round(bardata$Flux_Sum,0),")",sep=""), )+
    geom_vline(xintercept=0, color="Gray20", linewidth = 2)
  #ggsave("Figure4 - (b) Korea.png", units="in", width=6, height=6, dpi=1500)
}
for(plot in "Figure4-c"){ 
  #[Chl-a]
  # Korea
  bardata=read.csv("Data-Fig4.csv") %>% filter(Species=="Chl_a") %>% filter(Region=="Korean")
  ggplot()+
    geom_col_pattern(data=bardata, aes(x=Flux_Sum_percent, y=Type, pattern=Type, pattern_angle=Type),
                     fill="white",color="black", width=0.7, linewidth=2,
                     pattern_density=0.05, pattern_fill="Grey50", pattern_color="Grey50",
                     pattern_spacing=c(0.04,0.02,0.05, 1))+
    scale_y_discrete(limits=c("error", "dQdt_C", "dCdt_Q", "dCQdt"))+
    coord_cartesian(xlim = c(-100,50))+ 
    theme_void()+
    theme(legend.position = "none", axis.ticks.length = unit(0,"mm"),
          axis.text = element_blank(), axis.ticks = element_blank())+
    annotate("text", x = 34, y = bardata$Type, angle = -90, color= "black", fontface=2, size=10,
             label = c("F*", "C*", "Q*", "E*"))+
    annotate("text", x = 22, y = bardata$Type, angle = -90, color= "black", fontface=2, size=10,
             label = paste(round(bardata$Flux_Sum_percent,1),"%",sep=""))+
    annotate("text", x = 10, y = bardata$Type, angle = -90, color= "black", fontface=2, size=8,
             label = paste("(",round(bardata$Flux_Sum,1),")",sep=""), )+
    geom_vline(xintercept=0, color="Gray20", linewidth = 2)
  #ggsave("Figure4 - (c) Korea.png", units="in", width=6, height=6, dpi=1500)
}
for(plot in "Figure4-d"){ 
  # [TN]
  # West 
  bardata <- read.csv("Data-Fig4.csv") %>% filter(Species=="TN", Region=="West")
  ggplot()+
    geom_col_pattern(data=bardata, aes(x=Flux_Sum_percent, y=Type, pattern=Type, pattern_angle=Type),
                     fill="white",color="black", width=0.7, linewidth=2,
                     pattern_density=0.05, pattern_fill="Grey50", pattern_color="Grey50",
                     pattern_spacing=c(0.04,0.02,0.05, 1))+
    scale_y_discrete(limits=c("error", "dQdt_C", "dCdt_Q", "dCQdt"))+
    coord_cartesian(xlim = c(50,-100))+
    theme_void()+
    theme(legend.position = "none", axis.ticks.length = unit(0,"mm"),
          axis.text = element_blank(), axis.ticks = element_blank())+
    annotate("text", x = 25, y = bardata$Type, 
             label = paste(c(round(bardata$Flux_Sum_percent[1],1),
                             round(bardata$Flux_Sum_percent[2],1),
                             round(bardata$Flux_Sum_percent[3],0),
                             round(bardata$Flux_Sum_percent[4],1)),
                           "%",sep=""), 
             color= "black", fontface=2, size=15)+
    geom_vline(xintercept=0, color="black", linewidth = 2)
  
  #ggsave("Figure4 - (d) West.png", units="in", width=7, height=7, dpi=900)
  
  # South
  bardata <- read.csv("Data-Fig4.csv") %>% filter(Species=="TN", Region=="South")
  ggplot()+
    geom_col_pattern(data=bardata, aes(x=Flux_Sum_percent, y=Type, pattern=Type, pattern_angle=Type),
                     fill="white",color="black", width=0.7, linewidth=2,
                     pattern_density=0.05, pattern_fill="Grey50", pattern_color="Grey50",
                     pattern_spacing=c(0.04,0.02,0.05, 1))+
    scale_y_discrete(limits=c("error", "dQdt_C", "dCdt_Q", "dCQdt"))+
    coord_cartesian(xlim = c(80,-70))+ 
    theme_void()+
    theme(legend.position = "none", axis.ticks.length = unit(0,"mm"),
          axis.text = element_blank(), axis.ticks = element_blank())+
    annotate("text", x = 10, y = bardata$Type, angle=-90,
             label = paste(c(round(bardata$Flux_Sum_percent[1],0),
                             round(bardata$Flux_Sum_percent[2],1),
                             round(bardata$Flux_Sum_percent[3],0),
                             round(bardata$Flux_Sum_percent[4],1)),
                           "%",sep=""), 
             color= "black", fontface=2, size=15)+
    geom_vline(xintercept=0, color="black", linewidth = 2)
  
  #ggsave("Figure4 - (d) South.png", units="in", width=7, height=7, dpi=900)
  
  # East
  bardata <- read.csv("Data-Fig4.csv") %>% filter(Species=="TN", Region=="East")
  ggplot()+
    geom_col_pattern(data=bardata, aes(x=Flux_Sum_percent, y=Type, pattern=Type, pattern_angle=Type),
                     fill="white",color="black", width=0.7, linewidth=2,
                     pattern_density=0.05, pattern_fill="Grey50", pattern_color="Grey50",
                     pattern_spacing=c(0.04,0.02,0.05, 1))+
    scale_y_discrete(limits=c("error", "dQdt_C", "dCdt_Q", "dCQdt"))+
    coord_cartesian(xlim = c(-60,90))+ 
    theme_void()+
    theme(legend.position = "none", axis.ticks.length = unit(0,"mm"),
          axis.text = element_blank(), axis.ticks = element_blank())+
    annotate("text", x = bardata$Flux_Sum_percent + 25, y = bardata$Type, 
             label = paste(c(round(bardata$Flux_Sum_percent[1],0),
                             round(bardata$Flux_Sum_percent[2],1),
                             round(bardata$Flux_Sum_percent[3],1),
                             round(bardata$Flux_Sum_percent[4],1)),
                           "%",sep=""), 
             color= "black", fontface=2, size=15)+
    geom_vline(xintercept=0, color="black", linewidth = 2)
  
  #ggsave("Figure4 - (d) East.png", units="in", width=7, height=7, dpi=900)
}
for(plot in "Figure4-e"){ 
  # [TP]
  # West 
  bardata <- read.csv("Data-Fig4.csv") %>% filter(Species=="TP", Region=="West")
  ggplot()+
    geom_col_pattern(data=bardata, aes(x=Flux_Sum_percent, y=Type, pattern=Type, pattern_angle=Type),
                     fill="white",color="black", width=0.7, linewidth=2,
                     pattern_density=0.05, pattern_fill="Grey50", pattern_color="Grey50",
                     pattern_spacing=c(0.04,0.02,0.05, 1))+
    scale_y_discrete(limits=c("error", "dQdt_C", "dCdt_Q", "dCQdt"))+
    coord_cartesian(xlim = c(50,-100))+
    theme_void()+
    theme(legend.position = "none", axis.ticks.length = unit(0,"mm"),
          axis.text = element_blank(), axis.ticks = element_blank())+
    annotate("text", x = 25, y = bardata$Type, 
             label = paste(c(round(bardata$Flux_Sum_percent[1],0),
                             round(bardata$Flux_Sum_percent[2],1),
                             round(bardata$Flux_Sum_percent[3],1),
                             round(bardata$Flux_Sum_percent[4],1)),
                           "%",sep=""), 
             color= "black", fontface=2, size=15)+
    geom_vline(xintercept=0, color="black", linewidth = 2)
  
  #ggsave("Figure4 - (e) West.png", units="in", width=7, height=7, dpi=900)
  
  # South
  bardata <- read.csv("Data-Fig4.csv") %>% filter(Species=="TP", Region=="South")
  ggplot()+
    geom_col_pattern(data=bardata, aes(x=Flux_Sum_percent, y=Type, pattern=Type, pattern_angle=Type),
                     fill="white",color="black", width=0.7, linewidth=2,
                     pattern_density=0.05, pattern_fill="Grey50", pattern_color="Grey50",
                     pattern_spacing=c(0.04,0.02,0.05, 1))+
    scale_y_discrete(limits=c("error", "dQdt_C", "dCdt_Q", "dCQdt"))+
    coord_cartesian(xlim = c(80,-70))+ 
    theme_void()+
    theme(legend.position = "none", axis.ticks.length = unit(0,"mm"),
          axis.text = element_blank(), axis.ticks = element_blank())+
    annotate("text", x = 10, y = bardata$Type, angle=-90,
             label = paste(c(round(bardata$Flux_Sum_percent[1],0),
                             round(bardata$Flux_Sum_percent[2],1),
                             round(bardata$Flux_Sum_percent[3],1),
                             round(bardata$Flux_Sum_percent[4],1)),
                           "%",sep=""), 
             color= "black", fontface=2, size=15)+
    geom_vline(xintercept=0, color="black", linewidth = 2)
  
  #ggsave("Figure4 - (e) South.png", units="in", width=7, height=7, dpi=900)
  
  # East
  bardata <- read.csv("Data-Fig4.csv") %>% filter(Species=="TP", Region=="East")
  ggplot()+
    geom_col_pattern(data=bardata, aes(x=Flux_Sum_percent, y=Type, pattern=Type, pattern_angle=Type),
                     fill="white",color="black", width=0.7, linewidth=2,
                     pattern_density=0.05, pattern_fill="Grey50", pattern_color="Grey50",
                     pattern_spacing=c(0.04,0.02,0.05, 1))+
    scale_y_discrete(limits=c("error", "dQdt_C", "dCdt_Q", "dCQdt"))+
    coord_cartesian(xlim = c(-60,90))+ 
    theme_void()+
    theme(legend.position = "none", axis.ticks.length = unit(0,"mm"),
          axis.text = element_blank(), axis.ticks = element_blank())+
    annotate("text", x = bardata$Flux_Sum_percent + 25, y = bardata$Type, 
             label = paste(c(round(bardata$Flux_Sum_percent[1],0),
                             round(bardata$Flux_Sum_percent[2],1),
                             round(bardata$Flux_Sum_percent[3],1),
                             round(bardata$Flux_Sum_percent[4],1)),
                           "%",sep=""),
             color= "black", fontface=2, size=15)+
    geom_vline(xintercept=0, color="black", linewidth = 2)
  
  #ggsave("Figure4 - (e) East.png", units="in", width=7, height=7, dpi=900)
}
for(plot in "Figure4-f"){ 
  #[Chl-a]
  # West 
  bardata <- read.csv("Data-Fig4.csv") %>% filter(Species=="Chl_a", Region=="West")
  ggplot()+
    geom_col_pattern(data=bardata, aes(x=Flux_Sum_percent, y=Type, pattern=Type, pattern_angle=Type),
                     fill="white",color="black", width=0.7, linewidth=2,
                     pattern_density=0.05, pattern_fill="Grey50", pattern_color="Grey50",
                     pattern_spacing=c(0.04,0.02,0.05, 1))+
    scale_y_discrete(limits=c("error", "dQdt_C", "dCdt_Q", "dCQdt"))+
    coord_cartesian(xlim = c(60,-90)+2)+
    theme_void()+
    theme(legend.position = "none", axis.ticks.length = unit(0,"mm"),
          axis.text = element_blank(), axis.ticks = element_blank())+
    annotate("text", x = c(20,35,20,20)+5, y = bardata$Type, 
             label = paste(c(signif(bardata$Flux_Sum_percent,2)), "%",sep=""), 
             color= "black", fontface=2, size=15)+
    geom_vline(xintercept=0, color="black", linewidth = 2)
  
  #ggsave("Figure4 - (f) West.png", units="in", width=7, height=7, dpi=900)
  
  # South
  bardata <- read.csv("Data-Fig4.csv") %>% filter(Species=="Chl_a", Region=="South")
  ggplot()+
    geom_col_pattern(data=bardata, aes(x=Flux_Sum_percent, y=Type, pattern=Type, pattern_angle=Type),
                     fill="white",color="black", width=0.7, linewidth=2,
                     pattern_density=0.05, pattern_fill="Grey50", pattern_color="Grey50",
                     pattern_spacing=c(0.04,0.02,0.05, 1))+
    scale_y_discrete(limits=c("error", "dQdt_C", "dCdt_Q", "dCQdt"))+
    coord_cartesian(xlim = c(80,-70))+ 
    theme_void()+
    theme(legend.position = "none", axis.ticks.length = unit(0,"mm"),
          axis.text = element_blank(), axis.ticks = element_blank())+
    annotate("text", x = 10, y = bardata$Type, angle=-90,
             label = paste(signif(bardata$Flux_Sum_percent,2), "%",sep=""), 
             color= "black", fontface=2, size=15)+
    geom_vline(xintercept=0, color="black", linewidth = 2)
  
  #ggsave("Figure4 - (f) South.png", units="in", width=8, height=8, dpi=900)
  
  # East
  bardata <- read.csv("Data-Fig4.csv") %>% filter(Species=="Chl_a", Region=="East")
  ggplot()+
    geom_col_pattern(data=bardata, aes(x=Flux_Sum_percent, y=Type, pattern=Type, pattern_angle=Type),
                     fill="white",color="black", width=0.7, linewidth=2,
                     pattern_density=0.05, pattern_fill="Grey50", pattern_color="Grey50",
                     pattern_spacing=c(0.04,0.02,0.05, 1))+
    scale_y_discrete(limits=c("error", "dQdt_C", "dCdt_Q", "dCQdt"))+
    coord_cartesian(xlim = c(-60,90))+ 
    theme_void()+
    theme(legend.position = "none", axis.ticks.length = unit(0,"mm"),
          axis.text = element_blank(), axis.ticks = element_blank())+
    annotate("text", x = bardata$Flux_Sum_percent + 25, y = bardata$Type, 
             label = paste(round(bardata$Flux_Sum_percent,1),"%",sep=""), 
             color= "black", fontface=2, size=15)+
    geom_vline(xintercept=0, color="black", linewidth = 2)
  
  #ggsave("Figure4 - (f) East.png", units="in", width=7, height=7, dpi=900)
}

for(plot in "Figure5-a"){
  #[TN]
  dual_axis_ratio = 5
  ggplot()+
    geom_errorbar(data=read.csv("Data-Fig5.csv") %>% filter(Species=="TN"),
                  aes(x=Xaxis, ymin=min, ymax=max, color=Type), width=0, linewidth=1)+
    geom_point(data=read.csv("Data-Fig5.csv") %>% filter(Species=="TN"),
               aes(x=Xaxis, y=median, shape=Type, fill=Type, size=Type), color="black", stroke=1)+
    scale_shape_manual(values=c("River('02-11)"=21, "River('12-21)"=24,
                                "Marginal Sea('02-11)"=19, "Marginal Sea('12-21)"=17),
                       limits=c("River('02-11)", "River('12-21)",
                                "Marginal Sea('02-11)", "Marginal Sea('12-21)"),
                       guide=guide_legend(ncol=2))+
    scale_fill_manual(values=c("River('02-11)"="blue3", "River('12-21)"="blue3",
                               "Marginal Sea('02-11)"="black", "Marginal Sea('12-21)"="black"),
                      limits=c("River('02-11)", "River('12-21)",
                               "Marginal Sea('02-11)", "Marginal Sea('12-21)"),
                      guide=guide_legend(ncol=2))+
    scale_color_manual(values=c("River('02-11)"="blue3", "River('12-21)"="blue3",
                                "Marginal Sea('02-11)"="black", "Marginal Sea('12-21)"="black"),
                       limits=c("River('02-11)", "River('12-21)",
                                "Marginal Sea('02-11)", "Marginal Sea('12-21)"),
                       guide=guide_legend(ncol=2))+
    scale_size_manual(values=c("River('02-11)"=3, "River('12-21)"=3,
                               "Marginal Sea('02-11)"=3, "Marginal Sea('12-21)"=3),
                      limits=c("River('02-11)", "River('12-21)",
                               "Marginal Sea('02-11)", "Marginal Sea('12-21)"),
                      guide=guide_legend(ncol=2))+
    
    geom_vline(xintercept = c(3,9,15), linetype="longdash", linewidth=0.5)+
    geom_vline(xintercept = c(6,12), linetype="solid", linewidth=0.8)+
    
    scale_x_continuous(breaks = c(3,9,15), labels=c("West", "South", "East"))+
    scale_y_continuous(breaks = seq(0, 10, 4),
                       sec.axis = sec_axis(~./dual_axis_ratio, name=bquote("Marginal Sea"),
                                           breaks=seq(0, 2, 0.8)))+
    coord_cartesian(ylim = c(0, 9))+
    xlab(NULL) + ylab(bquote("River"))+ 
    ggtitle(bquote("(a) TN ["~"g/m"^3~"]")) +
    theme_bw() + theme(panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(),
                       panel.grid.major.y = element_line(colour="gray92", size=0.5),
                       panel.grid.minor.y = element_line(colour="gray92", size=0.5),
                       axis.text.x = element_text(margin=margin(t=5), size=20, color="black"),
                       axis.title.x = element_text(margin=margin(t=10), size=20, color="black"),
                       axis.text.y.left = element_text(margin=margin(l=8,r=3), size=20, color="blue3"),
                       axis.title.y.left = element_text(margin = margin(r=3), size=20, color="blue3"), 
                       axis.text.y.right = element_text(margin=margin(r=10,l=3), size=20, color="black"),
                       axis.title.y.right  = element_text(margin = margin(l=3), size=20, color="black"), 
                       plot.title = element_text(margin=margin(b=15), hjust = 0, size = 20),
                       axis.ticks.length.y = unit(1.5,"mm"),axis.ticks.length.x = unit(0,"mm"),
                       axis.ticks.y = element_line(linewidth = 1), 
                       legend.position = c(0.732,0.855),
                       legend.background = element_rect(color="black", fill="white", linewidth=1),
                       legend.title = element_blank(), legend.text = element_text(size=12, face="bold"),
                       legend.key = element_rect(color=NA, fill=NA, linewidth=1),
                       legend.key.size=unit(5,"mm"), legend.direction = "horizontal",
                       legend.margin = margin(b=3, r=2), 
                       panel.border = element_rect(color="black", fill=NA, linewidth=1),
                       plot.margin = unit(c(1, 0.2, 0.2, 0.2),"cm"))
  
  #ggsave("Figure 5 - (a).png", units="in", width=8, height=3, dpi=1200)
}
for(plot in "Figure5-b"){
  dual_axis_ratio = 5
  ggplot()+
    geom_errorbar(data=read.csv("Data-Fig5.csv") %>% filter(Species=="TP"),
                  aes(x=Xaxis, ymin=min*10, ymax=max*10, color=Type), width=0, linewidth=1)+
    geom_point(data=read.csv("Data-Fig5.csv") %>% filter(Species=="TP"),
               aes(x=Xaxis, y=median*10, shape=Type, fill=Type, size=Type), color="black", stroke=1)+
    scale_shape_manual(values=c("River('02-11)"=21, "River('12-21)"=24,
                                "Marginal Sea('02-11)"=19, "Marginal Sea('12-21)"=17),
                       limits=c("River('02-11)", "River('12-21)",
                                "Marginal Sea('02-11)", "Marginal Sea('12-21)"),
                       guide=guide_legend(ncol=2))+
    scale_fill_manual(values=c("River('02-11)"="blue3", "River('12-21)"="blue3",
                               "Marginal Sea('02-11)"="black", "Marginal Sea('12-21)"="black"),
                      limits=c("River('02-11)", "River('12-21)",
                               "Marginal Sea('02-11)", "Marginal Sea('12-21)"),
                      guide=guide_legend(ncol=2))+
    scale_color_manual(values=c("River('02-11)"="blue3", "River('12-21)"="blue3",
                                "Marginal Sea('02-11)"="black", "Marginal Sea('12-21)"="black"),
                       limits=c("River('02-11)", "River('12-21)",
                                "Marginal Sea('02-11)", "Marginal Sea('12-21)"),
                       guide=guide_legend(ncol=2))+
    scale_size_manual(values=c("River('02-11)"=3, "River('12-21)"=3,
                               "Marginal Sea('02-11)"=3, "Marginal Sea('12-21)"=3),
                      limits=c("River('02-11)", "River('12-21)",
                               "Marginal Sea('02-11)", "Marginal Sea('12-21)"),
                      guide=guide_legend(ncol=2))+
    
    geom_vline(xintercept = c(3,9,15), linetype="longdash", linewidth=0.5)+
    geom_vline(xintercept = c(6,12), linetype="solid", linewidth=0.8)+
    
    scale_x_continuous(breaks = c(3,9,15), labels=c("West", "South", "East"))+
    scale_y_continuous(breaks = seq(0, 10, 2),
                       sec.axis = sec_axis(~./dual_axis_ratio, name=bquote("Marginal Sea"),
                                           breaks=seq(0, 2, 0.4)))+
    coord_cartesian(ylim = c(0, 5))+
    xlab(NULL) + ylab(bquote("River"))+ 
    ggtitle(bquote("(b) TP [10"^-1~"g/m"^3~"]")) +
    theme_bw() + theme(panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(),
                       panel.grid.major.y = element_line(colour="gray92", size=0.5),
                       panel.grid.minor.y = element_line(colour="gray92", size=0.5),
                       axis.text.x = element_text(margin=margin(t=5), size=20, color="black"),
                       axis.title.x = element_text(margin=margin(t=10), size=20, color="black"),
                       axis.text.y.left = element_text(margin=margin(l=8,r=3), size=20, color="blue3"),
                       axis.title.y.left = element_text(margin = margin(r=3), size=20, color="blue3"), 
                       axis.text.y.right = element_text(margin=margin(r=10,l=3), size=20, color="black"),
                       axis.title.y.right  = element_text(margin = margin(l=3), size=20, color="black"), 
                       plot.title = element_text(margin=margin(b=15), face = "bold", hjust = 0, size = 20),
                       axis.ticks.length.y = unit(1.5,"mm"),axis.ticks.length.x = unit(0,"mm"),
                       axis.ticks.y = element_line(linewidth = 1), 
                       legend.position = "none",legend.background = element_blank(),
                       legend.title = element_blank(), legend.text = element_text(size=15, face="bold"),
                       legend.key = element_rect(color=NA, fill=NA, linewidth=1),
                       legend.key.size=unit(6,"mm"), legend.direction = "horizontal",
                       legend.margin = margin(b=5), 
                       panel.border = element_rect(color="black", fill=NA, linewidth=1),
                       plot.margin = unit(c(1, 0.2, 0.2, 0.2),"cm"))
  
  #ggsave("Figure 5 - (b).png", units="in", width=8, height=3, dpi=1200)
}
for(plot in "Figure5-c"){
  dual_axis_ratio = 5
  ggplot()+
    geom_errorbar(data=read.csv("Data-Fig5.csv") %>% filter(Species=="Chl-a"),
                  aes(x=Xaxis, ymin=min*100, ymax=max*100, color=Type), width=0, linewidth=1)+
    geom_point(data=read.csv("Data-Fig5.csv") %>% filter(Species=="Chl-a"),
               aes(x=Xaxis, y=median*100, shape=Type, fill=Type, size=Type), color="black", stroke=1)+
    scale_shape_manual(values=c("River('02-11)"=21, "River('12-21)"=24,
                                "Marginal Sea('02-11)"=19, "Marginal Sea('12-21)"=17),
                       limits=c("River('02-11)", "River('12-21)",
                                "Marginal Sea('02-11)", "Marginal Sea('12-21)"),
                       guide=guide_legend(ncol=2))+
    scale_fill_manual(values=c("River('02-11)"="blue3", "River('12-21)"="blue3",
                               "Marginal Sea('02-11)"="black", "Marginal Sea('12-21)"="black"),
                      limits=c("River('02-11)", "River('12-21)",
                               "Marginal Sea('02-11)", "Marginal Sea('12-21)"),
                      guide=guide_legend(ncol=2))+
    scale_color_manual(values=c("River('02-11)"="blue3", "River('12-21)"="blue3",
                                "Marginal Sea('02-11)"="black", "Marginal Sea('12-21)"="black"),
                       limits=c("River('02-11)", "River('12-21)",
                                "Marginal Sea('02-11)", "Marginal Sea('12-21)"),
                       guide=guide_legend(ncol=2))+
    scale_size_manual(values=c("River('02-11)"=3, "River('12-21)"=3,
                               "Marginal Sea('02-11)"=3, "Marginal Sea('12-21)"=3),
                      limits=c("River('02-11)", "River('12-21)",
                               "Marginal Sea('02-11)", "Marginal Sea('12-21)"),
                      guide=guide_legend(ncol=2))+
    
    geom_vline(xintercept = c(3,9,15), linetype="longdash", linewidth=0.5)+
    geom_vline(xintercept = c(6,12), linetype="solid", linewidth=0.8)+
    
    scale_x_continuous(breaks = c(3,9,15), labels=c("West", "South", "East"))+
    scale_y_continuous(breaks = seq(0, 10, 2),
                       sec.axis = sec_axis(~./dual_axis_ratio, name=bquote("Marginal Sea"),
                                           breaks=seq(0, 2, 0.4)))+
    coord_cartesian(ylim = c(0, 5.5))+
    xlab(NULL) + ylab(bquote("River"))+ 
    ggtitle(bquote("(c) Chl-a [10"^-2~"g/m"^3~"]")) +
    theme_bw() + theme(panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(),
                       panel.grid.major.y = element_line(colour="gray92", size=0.5),
                       panel.grid.minor.y = element_line(colour="gray92", size=0.5),
                       axis.text.x = element_text(margin=margin(t=5), size=20, color="black"),
                       axis.title.x = element_text(margin=margin(t=10), size=20, color="black"),
                       axis.text.y.left = element_text(margin=margin(l=8,r=3), size=20, color="blue3"),
                       axis.title.y.left = element_text(margin = margin(r=3), size=20, color="blue3"), 
                       axis.text.y.right = element_text(margin=margin(r=10,l=3), size=20, color="black"),
                       axis.title.y.right  = element_text(margin = margin(l=3), size=20, color="black"), 
                       plot.title = element_text(margin=margin(b=15), face = "bold", hjust = 0, size = 20),
                       axis.ticks.length.y = unit(1.5,"mm"),axis.ticks.length.x = unit(0,"mm"),
                       axis.ticks.y = element_line(linewidth = 1), 
                       legend.position = "none",legend.background = element_blank(),
                       legend.title = element_blank(), legend.text = element_text(size=15, face="bold"),
                       legend.key = element_rect(color=NA, fill=NA, linewidth=1),
                       legend.key.size=unit(6,"mm"), legend.direction = "horizontal",
                       legend.margin = margin(b=5), 
                       panel.border = element_rect(color="black", fill=NA, linewidth=1),
                       plot.margin = unit(c(1, 0.2, 0.2, 0.2),"cm"))
  
  #ggsave("Figure 5 - (c).png", units="in", width=8, height=3, dpi=1200)
}

for(plot in "Figure6-a"){
  ggplot()+
    geom_errorbar(data=read.csv("Data-Fig6(a).csv") %>% filter(Xaxis %in% sort(c(seq(1,62,3),seq(2,62,3)))), 
                  aes(x=Xaxis, ymin=min, ymax=max), width=0.8, color="black", linewidth=0.75)+
    geom_point(data=read.csv("Data-Fig6(a).csv") %>% filter(Xaxis %in% sort(c(seq(1,62,3),seq(2,62,3)))), 
               aes(x=Xaxis, y=median, shape=Decade), size=3, color="black")+
    
    geom_vline(xintercept = seq(3,62,3), linetype="dashed", linewidth=0.5)+
    
    scale_shape_manual('Decade Mean', labels = c('P'="2002-2011", 'R'="2012-2021"),
                       values = c('P'=16, 'R'=17))+
    scale_x_continuous(breaks = sort(seq(1.5, 62, 3)),
                       labels = paste0("W",sprintf("%02d", sort(seq(1,21,1)))))+
    scale_y_continuous(breaks = seq(0, 200, 50))+
    coord_cartesian(ylim = c(-5, 205), xlim = c(0,63), expand = FALSE)+
    ylab(NULL)+ xlab(NULL) + 
    ggtitle(bquote("(a) Western River Chl-a ["~"mg/m"^3~"]"))+
    theme_bw() + theme(panel.grid.major.x = element_line(colour="white"),
                       panel.grid.minor.x = element_line(colour="white"),
                       axis.text.x = element_text(margin=margin(t=5), size=15, color="black"),
                       axis.text.y.left = element_text(margin=margin(l=8,r=3), size=18, color="black"),
                       axis.title.y.left = element_text(margin = margin(r=3), size=18, color="black"), 
                       plot.title = element_text(margin=margin(b=10), face = "bold", hjust = 0, size = 20),
                       axis.ticks.length.y = unit(1.5,"mm"),axis.ticks.length.x = unit(0,"mm"),
                       axis.ticks.y = element_line(linewidth = 1), 
                       legend.position = "none",
                       panel.border = element_rect(color="black", fill=NA, linewidth=1))
  
  #ggsave("Figure 6 - (a).png", units="in", width=16, height=4, dpi=1200)
}
for(plot in "Figure6-b"){
  ggplot()+
    geom_errorbar(data=read.csv("Data-Fig6(b).csv") %>% filter(Xaxis %in% sort(c(seq(1,114,3),seq(2,114,3)))), 
                  aes(x=Xaxis, ymin=min, ymax=max), width=0.8, color="black", linewidth=0.75)+
    geom_point(data=read.csv("Data-Fig6(b).csv") %>% filter(Xaxis %in% sort(c(seq(1,114,3),seq(2,114,3)))), 
               aes(x=Xaxis, y=median, shape=Decade), size=3, color="black")+
    
    geom_vline(xintercept = seq(3,114,3), linetype="dashed", linewidth=0.5)+
    
    scale_shape_manual('Decade Mean', labels = c('P'="2002-2011", 'R'="2012-2021"),
                       values = c('P'=16, 'R'=17))+
    scale_x_continuous(breaks = sort(seq(1.5, 114, 3)),
                       labels = paste0("M",sprintf("%02d", sort(seq(1,38,1)))))+
    scale_y_continuous(breaks = seq(0, 20, 5))+
    coord_cartesian(ylim = c(-0.5, 21), xlim = c(0,115), expand = FALSE)+
    ylab(NULL)+ xlab(NULL) +
    ggtitle(bquote("(b) Western Marginal Sea Chl-a ["~"mg/m"^3~"]"))+
    theme_bw() + theme(panel.grid.major.x = element_line(colour="white"),
                       panel.grid.minor.x = element_line(colour="white"),
                       axis.text.x = element_text(margin=margin(t=5), size=13, color="black"),
                       axis.text.y.left = element_text(margin=margin(l=8,r=3), size=18, color="black"),
                       axis.title.y.left = element_text(margin = margin(r=3), size=18, color="black"), 
                       plot.title = element_text(margin=margin(b=10), hjust = 0, size = 20),
                       axis.ticks.length.y = unit(1.5,"mm"),axis.ticks.length.x = unit(0,"mm"),
                       axis.ticks.y = element_line(linewidth = 1), 
                       legend.position = "none",
                       panel.border = element_rect(color="black", fill=NA, linewidth=1))
  
  #ggsave("Figure 6 - (b).png", units="in", width=16, height=4, dpi=1200)
}